# Figure2_Long-term 12 month SPI values.R
# Coded by Prof. Alastair J. Potts
# 31 Jan 2025
##########################################
# Use groundhog package for reproducibility
use.groundhog <- TRUE

if (use.groundhog) {
  library(groundhog)  # Load groundhog package
  set.groundhog.folder("C:\\Groundhog_R\\R4.4.1_2024-05-14")  # Set groundhog folder location
  GroundhogDay <- '2024-05-14'  # Define GroundhogDay version
  groundhog.library("tidyverse",GroundhogDay)  # Load required packages using groundhog
  groundhog.library("lubridate",GroundhogDay)
  groundhog.library("drought",GroundhogDay)
  groundhog.library("SPEI",GroundhogDay)
  groundhog.library("patchwork",GroundhogDay)
} else {
  library("tidyverse")  # Load required libraries normally if not using groundhog
  library("lubridate")
  library("drought")
  library("SPEI")
  library("patchwork")
}

# Set working directory
setwd("D:\\Dropbox\\100_PROJECTS\\2024_ThicketHydraulics_SkeltonButtnerPotts")

# Read rainfall data, skipping the first 3 rows
read_csv("1_data/02_MurrayParkRainfall.csv", skip=3) -> dat

# Data transformation
# Convert monthly rainfall from inches to millimeters, remove unnecessary columns, and reshape data to long format
dat %>%
  mutate(across(Jan:Dec, function(fx) {fx * 25.4})) %>%  # Convert inches to mm
  select(-UNITS, -Tot, -CalculatedTotal, -Comment) %>%  # Remove unwanted columns
  pivot_longer(cols = Jan:Dec, names_to = "Month", values_to = "ppt") %>%  # Reshape data to long format
  filter(!is.na(ppt)) %>%  # Remove missing values
  as.data.frame() ->
  tutu

# Compute 12-month Standardized Precipitation Index (SPI)
spi_val <- spi(as.numeric(tutu[,"ppt"]), 12) %>%
  .$fitted %>% data.frame()

# Merge SPI values and plot time series
tutu %>%
  bind_cols(spi_val) %>%
  as_tibble() %>%
  rename(spi_val = ".") %>%
  mutate(ym = ym(paste(Year, Month)),
         ym = ceiling_date(ym, unit = "month") - 1) %>%
  ggplot(aes(ym, spi_val)) +
  geom_hline(yintercept = c(-3, -2, -1, 1, 2, 3), linetype = "dashed", colour = "darkgrey") +
  geom_line(size = 0.5) +
  theme_bw() +
  scale_y_continuous(breaks = seq(-3, 3)) +
  scale_x_date(date_breaks = "year", date_labels = "%Y",
               limits = c(ymd("1922-01-01"), ymd("2022-12-31")),
               expand = c(0, 0)) +
  ylab("12-month SPI") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
        legend.position = "none") +
  xlab("") +
  geom_hline(yintercept = 0)

# Generate SPI time series for different time frames and save plots

# SPI for 1984-2023
tutu %>%
  bind_cols(spi_val) %>%
  as_tibble() %>%
  rename(spi_val = ".") %>%
  mutate(ym = ym(paste(Year, Month)),
         ym = ceiling_date(ym, unit = "month") - 1) %>%
  ggplot(aes(ym, spi_val)) +
  geom_hline(yintercept = c(-3, -2, -1, 1, 2, 3), linetype = "dashed", colour = "darkgrey") +
  geom_line(size = 0.5) +
  theme_bw() +
  scale_y_continuous(breaks = seq(-3, 3)) +
  scale_x_date(date_breaks = "year", date_labels = "%Y",
               limits = c(ymd("1984-01-01"), ymd("2022-12-31")),
               expand = c(0, 0)) +
  ylab("12-month SPI") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
        legend.position = "none") +
  xlab("") +
  geom_hline(yintercept = 0) ->
  MP_84

# SPI for 1922-2023 with a shaded area for pre-1984 period
tutu %>%
  bind_cols(spi_val) %>%
  as_tibble() %>%
  rename(spi_val = ".") %>%
  mutate(ym = ym(paste(Year, Month)),
         ym = ceiling_date(ym, unit = "month") - 1) %>%
  ggplot(aes(ym, spi_val)) +
  geom_hline(yintercept = c(-3, -2, -1, 1, 2, 3), linetype = "dashed", colour = "darkgrey") +
  annotate("rect", xmin = ymd("1922-01-01"), xmax = ymd("1984-01-01"),
           ymin = -3, ymax = 3,
           alpha = .5, fill = "grey") +
  geom_line(size = 0.5) +
  theme_bw() +
  scale_y_continuous(breaks = seq(-3, 3), expand = c(0, 0)) +
  scale_x_date(breaks = seq(from = ymd("1925-01-01"), to = ymd("2022-12-31"), by = "5 years"),
               date_labels = "%Y",
               limits = c(ymd("1922-01-01"), ymd("2022-12-31")),
               expand = c(0, 0)) +
  ylab("12-month SPI") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
        legend.position = "none") +
  xlab("") +
  geom_hline(yintercept = 0) ->
  MP_22

# Read CHIRPS dataset

read.csv("1_data/03_CHIRPS.csv",skip=3) -> chirps_spi

# SPI plot for CHIRPS dataset
chirps_spi %>%
  mutate(ym = ym(paste(Year, Month)),
         ym = ceiling_date(ym, unit = "month") - 1) %>%
  ggplot(aes(ym, spi_val)) +
  geom_hline(yintercept = c(-3, -2, -1, 1, 2, 3), linetype = "dashed", colour = "darkgrey") +
  geom_line(size = 0.5) +
  theme_bw() +
  scale_y_continuous(breaks = seq(-3, 3)) +
  scale_x_date(date_breaks = "year", date_labels = "%Y",
               limits = c(ymd("1984-01-01"), ymd("2022-12-31")),
               expand = c(0, 0)) +
  ylab("12-month SPI") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
        legend.position = "bottom") +
  xlab("") +
  geom_hline(yintercept = 0) ->
  fig_kab

#########################################################################
# Create the compound figure. 
(MP_22+
   ggtitle('A) Russel Park: 1922-2023'))/
  (MP_84+theme(axis.title.x=element_blank(),
               axis.text.x=element_blank(),
               axis.ticks.x=element_blank())+
     ggtitle('B) Russel Park: 1984-2023'))/
  (fig_kab+ggtitle('C) Kaboega: 1984-2023 (CHIRPS)'))

#########################################################################

ggsave("results/Figure2.jpg",units="cm",height=15,width=15,dpi = 800)
ggsave("results/Figure2.tif",units="cm",height=15,width=15,dpi = 800)
ggsave("results/Figure2.pdf",units="cm",height=15,width=15)

